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Abstract. We study the time evolution of two protoplanets still embedded in a protoplanetary disk. The results 
of two different numerical approaches are presented and compared. In the hrst approach, the motion of the disk 
material is computed with viscous hydrodynamical simulations, and the planetary motion is determined by N- 
body calculations including exactly the gravitational forces exerted by the disk material. In the second approach, 
only the N-body integration is performed but with additional dissipative forces included such as to mimic the 
effect of the disk torques acting on the disk. This type of modeling is much faster than the full hydrodynamical 
simulations, and gives comparative results provided that parameters are adjusted properly. 

Resonant capture of the planets is seen in both approaches, where the order of the resonance depends on the 
properties of the disk and the planets. Resonant capture leads to a rise in the eccentricity and to an alignment 
of the spatial orientation of orbits. The numerical results are compared with the observed planetary systems in 
mean motion resonance (Gl 867, HD 82943, and 55 Cnc). We find that the forcing together of two planets by their 
parent disk produces resonant configurations similar to those observed, but that eccentricity damping greater 
than that obtained in our hydrodynamic simulations is required to match the GJ 876 observations. 
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1. Introduction 

Since their first discovery in 1995, the number of detected 
extrasolar planets orbiting solar-type stars has risen dur- 
ing recent years to more than 100 (for an up-to-date list 
see e.g. http://www.obspm.fr/encycl/encycl.html by 
J. Schneider). Among these, there are currently 11 systems 
with two or more planets; a summary of their properties 
has been given recently by Marcy et al. (2003). With fur- 
ther observations to come, the fraction of systems with 
multiple planets will almost certainly increase, as many 
of the systems exhibit long-term trends in their radial ve- 
locity, suggesting an additional outer planet. Among the 
known multiple-planet extrasolar systems there are now 
three confirmed cases, namely Gl 876 (Marcy et al. 2001), 
HD 82943 (the Coralie Planet Search Programme, ESO 
Press Release 07/01), and 55 Cnc (Marcy et al. 2002) 
where the planets orbit their central star in a low-order 
mean motion resonance such that the orbital periods have 
nearly exactly the ratio 2:1 or 3:1. The parameters of these 
planetary systems are displayed in Table 1 below. The pos- 
sibility of a 2:1 resonance in HD 160691 has also been dis- 
cussed recently by Bois et al. (2003), although the orbital 
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periods are too long to definitely confirm this. Overall, 
these numbers imply that at least one-fourth of multiple- 
planetary systems contain planets in resonance, a fraction 
which is even higher if secular resonances, such as those 
observed in the v And system, (Butler et al. 1999) are also 
considered. 

The formation of resonant planetary systems can be 
understood by considering the joint evolution of proto- 
planets together with the protoplanetary disk from which 
they formed. Using local linear analysis, it has been shown 
that the gravitational interaction of a single protoplanet 
with its disk leads to torques resulting in a change of 
the semi-major axis (migration) of the planet (Goldreich 
& Tremaine 1980; Lin & Papaloizou 1986; Ward 1997; 
Tanaka et al. 2002). Additionally, as a result of angu- 
lar momentum transfer between the viscous disk and the 
planet, planetary masses of around one Jupiter mass can 
open gaps in the surrounding disk (Lin & Papaloizou 
1980, 1993). Fully non-linear hydrodynamical calculations 
for Jupiter-sized planets (Kley 1999; Bryden et al. 1999; 
Lubow et al. 1999; Nelson et al. 2000; D'Angelo et al. 2002) 
confirmed this expectation and clearly showed that disk- 
planet interaction leads to: i) excitation of spiral shock 
waves in the disk, whose tightness depends on the sound- 
speed in the disk, ii) formation of an annular gap, whose 
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width is determined by the balance between gap-opening 
tidal torques and gap-closing viscous plus pressure forces, 
Hi) inward migration on a time scale of 10 5 yrs for typical 
disk parameters, in particular disk masses corresponding 
to that of the minimum mass solar nebula, iv) possible 
mass growth after gap formation up to about 10 Mj up 
when finally the gravitational torques overwhelm the dif- 
fusive tendencies of the gas, and v) a prograde rotation 
of the planet. New three-dimensional computations with 
high resolution resolve the flow structure in the vicinity of 
the planet, and allow for more accurate estimates of the 
mass accretion and migration rates (D'Angelo et al. 2003; 
Bate et al. 2003). 

These hydrodynamic simulations with single planets 
have been extended to models which contain multiple 
planets. It has been shown (Kley 2000; Bryden et al. 2000; 
Snellgrove et al. 2001; Nelson & Papaloizou 2002) that 
during the early evolution, when the planets are still em- 
bedded in the disk, different migration speeds may lead 
to an approach of neighboring planets and eventually to 
resonant capture. More specifically, the evolution of plan- 
etary systems into a 2:1 resonant configuration was seen 
in the calculations of Kley (2000) prior to the discovery of 
any such systems. 

In addition to hydrodynamic disk-planet simulations, 
many authors have analyzed the evolution of multiple- 
planet systems with N-body methods. Each of the known 
resonant systems have been considered in detail. Ji et al. 
(2002) and Lee & Peale (2002a) have modeled the evo- 
lution of 2:1 resonant system GJ 876, while the 3:1 sys- 
tem 55 Cnc has been analyzed by Ji et al. (2003b) and 
Lee & Peale (2002b), and the 2:1 system HD 82943 by 
Gozdziewski & Maciejewski (2001) and Ji et al. (2003a). 
Based on orbit integrations, these papers confirm that 
the planets in these systems are in resonance with each 
other. The dynamics and stability of resonant planetary 
systems in general has been recently studied by Beauge 
et al. (2002). 

Here we present new numerical calculations treating 
the evolution of two planets still embedded in a proto- 
planetary disk. We use both hydrodynamical simulations 
and simplified N-body integrations to follow the evolution 
of the system. In the first approach, the disk is evolved by 
solving the full time-dependent Navicr-Stokcs equations 
simultaneously with the evolution of the planets. Here, 
the motion of the planets is determined by the gravita- 
tional action of both planets, the star, and the disk. In 
the latter approach, we take a simplified approximation 
and perform 3-body (star plus two planets) calculations 
augmented by additional (damping) forces which approx- 
imately account for the gravitational influence of the disk 
(e.g. Lee & Peale 2002a). Using both approaches, allows 
a direct comparison of the alternative methods, and does 
enable us to determine the damping parameters required 
for the simpler (and much faster) second type of approach. 



Table 1. The orbital parameters of the three systems known to 
contain a mean motion resonance. P denotes the orbital period, 
Msini the mass of the planets, a the semi-major axis, e the 
eccentricity, vj the angle of periastron, and M* the mass of the 
central star. It should be noted that the orbital elements for 
shorter period planets undergo secular time variations. Thus 
in principle one should always state the epoch corresponding 
to these osculating elements (see e.g. Laughlin & Chambers 
2001). 



Name 


P 


M sin i 


a 


e 


vj 


M* 




[d] 


[M Jup ] 


[AU] 




[deg] 


[M Q ] 


Gl 876 


(2:1) 










0.32 


c 


30.1 


0.56 


0.13 


0.24 


159 




b 


61.02 


1.89 


0.21 


0.04 


163 




HD 82943 


(2:1) 










1.05 


b 


221.6 


0.88 


0.73 


0.54 


138 




c 


444.6 


1.63 


1.16 


0.41 


96 




55 Cnc 


(3:1) 










0.95 


b 


14.65 


0.84 


0.11 


0.02 


99 




c 


44.26 


0.21 


0.24 


0.34 


61 




d 


5360 


4.05 


5.9 


0.16 


201 





2. The Observations 

The basic orbital parameters of the three known systems 
in mean motion resonance are presented in Table 1. The 
orbital parameters for Gl 876 are taken from the dynami- 
cal fit of Laughlin & Chambers (2001), and for HD 82943 
from Gozdziewski & Maciejewski (2001). Due to the un- 
certainty in the inclinations of the systems, Msini, rather 
than the exact mass of each planet, is listed. By including 
the mutual perturbations of the planets into their fit of 
Gl 876, Laughlin & Chambers (2001), however, are able 
to constraint that system's inclination to ~ 30° — 50°. 

Two of the systems, Gl 876 and HD 82943, arc in a 
nearly exact 2:1 resonance. We note that in both cases 
the outer planet is more massive, in one case by a factor 
of about two (HD 82943) and in the other by more than 
three (Gl 876). The eccentricity of the inner (less massive) 
planet is larger than that of the outer one in both systems. 
For the system Gl 876 the alignment of the orbits is such 
that the two periastrae are pointing in nearly the same 
direction. For the system HD 82943 these data have not 
been clearly identified, due to the much longer orbital pe- 
riods, but they do not seem to be very different from each 
other. The third system, 55 Cnc, is actually a triple sys- 
tem. Here the inner two planets orbit the star very closely 
and are in a 3:1 resonance, while the third, most massive 
planet orbits at a distance of several AU. 

3. The Models 

Our goal is to investigate the evolution of protoplanets 
still embedded in their disk. As outlined in the intro- 
duction, we employ two different methods which comple- 
ment each other. First, a fully time-dependent hydrody- 
namical model for the joint evolution of the planets and 
disk is presented. Because the evolutionary time scale may 
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cover several thousands of orbits, the fully hydrodynami- 
cal computations (of disk and planets) can require millions 
of time-steps, which translates into an effective computa- 
tional time of up to several weeks. 

Since our main interest is the orbital evolution of the 
planets and not so much the hydrodynamics of the disk, 
we also perform 3-body orbit integrations which do not 
explicitly follow the disk's evolution. Through a direct 
comparison with the hydrodynamical models, it is then 
possible to infer the effective damping forces to include 
within this faster calculation. 

3.1. Hydrodynamical Model 

The first set of coupled hydrodynamical-N-body models 
presented in this paper are calculated in the same manner 
as the models described previously in Kley (1998, 1999) 
for single planets and in Kley (2000) for multiple planets. 
The reader is referred to those papers for details on the 
computational aspects of the simulations. Other similar 
models, following explicitly the motion of single and mul- 
tiple planets in disks, have been presented by Nelson et al. 
(2000), Bryden et al. (2000), and Snellgrove et al. (2001). 

3.1.1. Equations 

For reference, we summarize the basic equations for this 
problem. We use two-dimensional cylindrical coordinates 
(r, (p) , where r is the radial coordinate and ip is the az- 
imuthal angle. Thus, we consider an infinitesimally thin 
disk located at z = 0, with a velocity field u = {u r ,u v ). 
The origin of the coordinate system is at the position of 
the star. In the following we will use the symbol v = u r for 
the radial velocity and lo = u v /r for the angular velocity 
of the flow. As there is no preferred rotational frame, we 
work in a non-rotating reference system. Then the equa- 
tions of motion are 



BY 

- + V.( S u)=0, 



d(Y,v) 



dt 

d(Er 2 uj) 
dt 



V • (Suu) = £ rcu 2 



dp 
dr 



+ V • (£r 2 um) 



dp 

dip 



dr 



(1) 



(2) 



Here £ denotes the surface density, p the two-dimensional 
pressure, and f r , f v denote the components of the viscous 
forces, given explicitly in Kley (1999). The gravitational 
potential <I> generated by the protostar with mass M* and 
the two planets having mass mi and m 2 is given by 



$ = - 



r-r*| [( r _ ri )2 + s 2]V2 



Gm 2 



[(r-r 2 )2 + *|] 1/2 



(4) 



where G is the gravitational constant and r„, ri, and r 2 
are the position vectors for the star and the inner/outer 
planet, respectively. The quantities s\ and s 2 are smooth- 
ing lengths that are set to 1/5 of each planet's Roche 



radius. This softening of the potential removes any local 
fluctuations that might result as the planets move through 
the computational grid. 

The motion of the star and the planets is determined 
firstly by their mutual gravitational interaction and sec- 
ondly by the gravitational forces exerted upon them by the 
disk. The acceleration of the star a* is given for example 

by 



a, = Gm\ 



ri 



- r . |3 



Gm 2 



r 2 - r* 



r 2 



-r.|3 



+ G 



£(r, ip) 



■ dA (5) 



Disk 



where the integration covers the entire disk surface. The 
accelerations of the planets are found similarly. We work 
here in an accelerated coordinate frame where the origin 
is located in the center of the (moving) star. Thus, in 
addition to force due to the gravitational potential (Eq. 4) , 
the disk and planets also feel an acceleration —a*. 

3.1.2. Initial and boundary conditions 

The initial hydrodynamic structure of the disk, which ex- 



tends radially from 



to r max , is axisymmetric with 



respect to the location of the star, and the surface density 
scales as S(r) = S r -1 / 2 , with superimposed initial gaps 
(Kley 2000). The initial velocity is pure Keplerian rotation 
(v r = 0, v v — GM^/r 1 / 2 ). We assume a fixed temperature 
law with T(r) oc r _1 which follows from the assumed con- 
stant vertical height H/r. The kinematic viscosity v is 
typically parameterized by an a-description v = ac s H, 
with the sound speed c s — Hv v /r, only one model (X) 
has a constant kinematic viscosity. 

The radial outer boundary is closed, i.e. n r (r mm ) = 0. 
At the inner radial boundary outflow boundary condi- 
tions are applied; matter may flow out, but none is al- 
lowed to enter. This procedure mimics the accretion pro- 
cess onto the star. The density gradient is set to zero at 
r min and r max , while the angular velocity there is fixed to 
be Keplerian. In the azimuthal direction, periodic bound- 
ary conditions for all variables are imposed. 



(3) 3.1.3. Model parameters 



We present several models that are listed in Table 2, for 
the complete model parameters see below. In all cases, 
the planets are allowed to migrate (change their semi- 
major axes) through the disk in accordance with the 
gravitational torques exerted on them. During the evolu- 
tion, material is removed from the centers of the planets' 
Roche-lobes and is assumed to have been accreted onto 
each planet, for the detailed procedure see Kley (1999). 
However, in spite of this assumed accretion process, in 
most models this mass is not added to the dynamical mass 
of the planets. Hence, the disk mass is slowly depleted 
while the planets' masses are held fixed at the initial val- 
ues. This assumption of constant planet mass throughout 
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Table 2. Planetary and disk parameters of the hydrodynamic 
models. The masses of the planets are given in Jupiter masses 
(M Jup = 1CP 3 M ). Index 1 refers to the inner and index 2 to 
the outer planet. For the last model X, we allow for variable 
planet masses which both start with 1 Mj up and then grow 
during the computation. For the viscosity, the corresponding 
value of a is given except for the last case (model X) in which 
we use a constant kinematic viscosity (corresponding to a = 
0.04 at a radial location of 5.2 AU). The relative disk scale 
height H/r is given in the last column. For the disk masses 
and computational domains, see section 3.1. 



Model 


mi 


7712 


Viscosity 


H/r 




[M Jup ] 


[Mj up ] 


a 




A 


5 


3 




0.10 


B 


5 


3 


2 • 10 -3 


0.10 


C 


3 


5 


10~' z 


0.10 


Dl, D2 


3 


5 


3 • 10~ 3 


0.10 


E 


3 


5 


10~ 3 


0.10 


F 


3 


5 


10~ 2 


0.075 


G 


3 


5 


10 -2 


0.05 


X 


1(V) 


1(V) 


v = const. 


0.05 



the computation is well justified, as the migration rate de- 
pends only weakly on the mass of the planet (Nelson et al. 
2000). Only in the last model X are the masses of the plan- 
ets allowed to grow during the computation, which may 
test this assumption. 

In all the models A-G, the disk extends radially from 
r min = 1 to r max = 30 AU. The disk mass within this 
radial range is about 0.04 M Q , and the planets are al- 
ways placed initially at 4 and 10 AU, respectively. In the 
first two models (A & B), the inner planet is more massive 
than the outer one, while in models C through G the inner 
planet is less massive. A range of values for the viscosity 
and disk thickness H/r are considered, as listed in Tab. 2. 
The values of a = 0.01 and H/r = 0.1 may be on the high 
side for protoplanetary disks but allow for a sufficiently 
rapid evolution of the system to identify clearly the gov- 
erning physical effects. The influence of these parameters 
is studied by comparison between different models. The 
parameters for models Dl and D2 are identical, but in the 
second case the density has been perturbed randomly by 
1%. 

Model X differs in several respects from the other mod- 
els. The kinematic viscosity v is constant and equivalent 
to a = 0.004 at 5.2 AU. Here, both embedded planets 
each have initial masses of 1 Mj up and are placed initially 
at 1 and 2 AU. This model is a continuation of the one 
presented previously in Kley (2000). The radial extent of 
the disk for this model was r m i„ = 1.3 to r max = 20.8 
AU, which contained a disk mass of 0.01M Q initially. The 
initial setup makes the surface density S(r) the same for 
all models. 

3.2. Damped N-body Model 



of orbits must be calculated. Hence, we also perform sim- 
pler N-body calculations of planetary systems with two 
planets orbiting a single star. We consider only coplanar 
systems, where the planetary orbits and the equatorial 
plane of the disk arc all aligned. The orbit of a single 
planet around a star is an ellipse and is determined by 
three orbital elements: the semi-major axis a, the eccen- 
tricity e, and the direction of periastron w. The actual 
location of the planet within the orbit can be obtained 
from the elapsed time since last periastron. In the case 
of a planetary system with more than one planet, due to 
the mutual gravitational perturbations the ellipses are no 
longer fixed in space. In the case where the masses of the 
planets arc much smaller that the stellar mass, at each 
epoch wc can fit an instantaneous ellipse to the orbit of 
each planet and obtain the corresponding osculating ele- 
ments of the orbit. These are calculated for each planet 
individually, considering only one planet at a time. 

The gravitational influence of the surrounding disk is 
modeled here through prescribed (damping) forces. We as- 
sume that these forces act on the momentary semi-major 
axis and eccentricity of the planets through explicitly spec- 
ified relations for a(t) and e(t), which vary with time. The 
changes d and e caused by the damping effects of the disk 
can be translated into additional forces changing directly 
the position x and velocity u of the planets. Our imple- 
mentation follows Lee & Peale (2002a), where explicit ex- 
pressions for these damping terms are given. As a test, we 
recalculated their model for GJ 876 and confirmed their 
results. 

Motivated by the basic idea of two planets orbiting 
inside a disk's cavity (see Fig. 1), we damp a and e for the 
outer planet only. We adopt a logarithmic time derivative 
of a of the form 

- = 77T with r(t) = r + fit, (6) 

a T{t) 

where r(t) denotes the damping time, and 3 a dimension- 
less positive 'stretching' constant. By making the ansatz 
Eq. (6) we tried to make the damping as simple as possi- 
ble with only two parameters to fit. In practice we found 
that the damping time r could not be chosen as a fixed 
constant, such that we assume a simple linear time depen- 
dence. Equation (6) can be integrated to yield 

a(t)=a Q -ll+ , (7) 

where o denotes the starting value of a at the initial time 
t = 0. 

The eccentricity damping is set to a fixed multiple of 
the semi-major axis damping 

B -=K a -, (8) 

e a 



As pointed out, the full hydrodynamical evolution is com- 
putationally very time consuming since tens of thousands 



where if is a constant, typically larger than unity (see also 
Lee & Peale 2002a). 
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Fig. 1. Overview of the density distribution of model C after 
1500 orbital periods of the inner planet. Higher density regions 
are brighter and lower ones are darker. The star lies at the cen- 
ter of the white inner region bounded by r m in = 1 AU. The 
location of the two planets is indicated by the white dots, and 
their Roche-lobes are also drawn. Clearly seen are the irregu- 
lar spiral wakes generated by the planets. Regular intertwined 
spiral arms are seen only outside of the second planet. 

The time scale To, the 'stretching' factor (3, and K 
are adjusted in order to match the results of the full hy- 
drodynamic calculations. Results of the two methods are 
compared in Section 4.2 below. 

4. Results 

The basic evolutionary sequence of two planets evolving 
simultaneously with a hydrodynamic disk has been cal- 
culated and described by Kley (2000) and Bryden et al. 
(2000). The model X was considered in Kley (2000), where 
it was found that both planets evolve into a 2:1 resonant 
configuration. Model X and model A are discussed fur- 
ther in a recent conference paper Kley (2003). Here, these 
models are listed essentially for completeness, while our 
main interest in this paper focuses on models C-G. 

To analyze the system dynamics in the presence of a 
disk, we monitor the evolution of the orbital elements a, e, 
and w of both planets throughout the simulations. In the 
case of a resonance, the orbits are coupled dynamically. 
For coplanar systems that are in a mean-motion commen- 
surability (p + q) :p, it suffices to use two resonant angles 
to describe the evolution. Here we consider the angular 
difference of the apsidal lines 

Aro = U72 — Wi (9) 



and the combination 

©i = (p + <l)^2 -pAi - qwi, (10) 

where \ are the mean longitudes of the inner (i = 1) and 
outer (i = 2) planets. The two resonant angles Aw and 
0i have also been used recently by Beauge et al. (2002) to 
study the possible stable solutions of resonant planetary 
systems. In Eq. (10), we have p = 1, q = 1 for a 2:1 
resonance, while for a 3:1 resonance p = 1, q = 2. 

Before we discuss details of resonant planetary evolu- 
tion we briefly summarize the main properties of model 
C, which serves as our standard reference model. 



4.1.1. Overview 

At the start of the simulations both planets are placed 
into an axisymmetric disk, where the density is initialized 
with partially opened gaps superimposed on an otherwise 
smooth radial density profile. Upon starting the evolution 
the two main effects arc: 

a) Because of the accretion of gas onto the two planets the 
radial region in between them is depleted in mass and 
finally cleared. This clearing time depends on the mass 
of the planets and on the viscosity and temperature of 
the disk. The individual gaps of higher mass planets 
are deeper, which lengthens the overall clearing tim- 
scale. Higher viscosity and temperature tend to 'push' 
material towards the planets and shorten the clearing 
time. Additionally, the disturbance by the two planets 
creates a strongly time dependent flow with two mu- 
tually interacting spiral arms which also pushes mat- 
ter towards the planets. The snapshot in Fig. 1 shows 
clearly this effect, which again reduces the clearing 
time. Thus, the high viscosity (a — 0.01) and tem- 
perature {H/R — 0.1) of the standard model (C) still 
allow for gap clearing in spite of the large mass of the 
planets. After about 5000 yrs, only 2% of the initial 
material between the planets remains. For model X, 
with lower mass planets, this phase takes only a few 
hundred orbital periods (Kley 2000). 
Concurrently with the central ring depletion, the re- 
gion interior to the inner planet loses material cither 
due to accretion onto the central star (the inner bound- 
ary is open to outflow) or by accretion onto the planet. 
As with the intermediate ring, the timescale for clear- 
ing this inner region again depends on the physical 
parameter of the system. Thus, after an initial tran- 
sient phase we typically expect the configuration of 
two planets orbiting within an inner cavity of the disk, 
as seen in Fig. 1 (see also Kley 2000). 

b) After initialization, the planets quickly (within a few 
orbital periods) excite non-axisymmetric disturbances, 
viz. the spiral waves, in the disk. In contrast to the 
single planet case these are not stationary in time, be- 
cause there is no preferred rotating frame. The grav- 
itational torques exerted on the two planets by those 
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non-axisymmetric density perturbations induce a mi- 
gration process for the planets. 
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Fig. 2. The semi-major axis (a), eccentricity (e) and position 
angle of the orbital periastron (pj) for the two planets ver- 
sus time for Model C. In this example, the planets have fixed 
masses of 3 and 5 Mj up , and are placed initially at 4 and 
10 AU, respectively. The inner planet is denoted by the black 
line, the outer by the light gray line. The dotted reference line 
(labeled 3:1), indicates the location of the 3:1 resonance with 
respect to the inner planet. 



Now, the planets' relative positions within the cavity have 
a distinct influence on their subsequent evolution. As a 
consequence of the clearing process, the inner planet is no 
longer surrounded by any disk material and thus cannot 
grow any further in mass. In addition, it cannot migrate 
anymore, because there is no torque-exciting material left 
in its vicinity. All the material of the outer disk is still 
available, on the other hand, to exert negative (Lindblad) 
torques on the outer planet. Hence, in the initial phase of 
the computations we observe an inwardly migrating outer 
planet and a stalled inner planet with a constant semi- 
major axis (see the first 5000 yrs in the top panel of Fig. 2). 

During the inward migration process the eccentricity of 
the outer planet remains small. As can be seen from Fig. 3 
there is always a sufficient amount of matter in the imme- 
diate vicinity (co-orbital region) of the planet to ensure 
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Fig. 3. The azimuthally averaged surface density (in dimen- 
sionless units) for model C at 40,000 yrs. The location of the 
semi-major axis of the planet is indicated by the vertical dashed 
line at r = 6.63, and the size of its Roche-radius by the solid 
circle. The left and right arrows indicate the locations of the 
inner and outer 2:1 Lindblad resonances, respectively. 



damping of the eccentricity. For a 5 Jupiter mass planet 
at r = 6.63AU the Roche-lobe size is 0.78 AU, which is 
indicated by the radius of the drawn circle. The arrows 
denote the position of the inner and outer 2:1 Lindblad 
resonances. 

The decrease in separation between the planets in- 
creases their gravitational interactions. Once the ratio of 
the planets' orbital periods has reached a ratio of two inte- 
gers, i.e. they are close to a mean motion resonance, reso- 
nant capture of the inner planet by the outer one may en- 
sue. Whether or not this does actually happen depends on 
the physical conditions in the disk (e.g. viscosity) and the 
orbital parameters of the planets. If the migration speed 
is too large, for example, there may not be enough time 
to excite the resonance, and the outer planet will con- 
tinue migrating inward (e.g. Haghighipour 1999). Also, if 
the initial eccentricities are too small, then there may be 
no capture, particularly for second-order resonances such 
as the 3:1 resonance (see e.g. Murray & Dermott 1999). 
For more details on capture probability see section (4.3.3) 
below. 



4.1.2. 3:1 Resonance: Model C 

The typical time evolution of the semi-major axis (a), ec- 
centricity (e) and direction of the periastron (w) are dis- 
played in Fig. 2 for the standard model C. The planets 
were initialized with zero eccentricities at distances of 4 
and 10 AU in a disk with partially cleared gaps. 
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Fig. 4. Evolution of Avj and Oi in model C. Left: The difference in the direction of the periastrae, Aza = zu2 — zui (in degrees) 
for both planets versus time (in thousand years). Right: The resonant angle for the 3:1 resonance Oi = 3A2 — Ai — 2roi (in 
degrees) as a function of time. 



In the beginning, after the inner gap has completely 
cleared, only the outer planet migrates inward, and the 
eccentricities of both planets remain relatively small, 
(< 0.02). After about 5000 yrs the outer planet has 
reached a semi-major axis with an orbital period three 
times that of the inner planet. The periodic gravitational 
forcing leads to the capture of the inner planet into a 3:1 
resonance with the outer one. This is indicated by the dot- 
ted reference line (labeled 3:1) in the top panel of Fig. (2), 
which marks the location of the 3:1 resonance with respect 
to the inner planet. 

We summarize the following important features of the 
evolution after resonant capture: 

a) In the course of the subsequent evolution, the outer 
planet, which is still driven inward by the outer 
disk material, forces the inner planet to also mi- 
grate inwards. Both planets migrate inward simulta- 
neously, always retaining their resonant configuration. 
Consequently, the migration speed of the outer planet 
slows down, and their radial separation declines. 

b) Upon resonant capture the eccentricities of both plan- 
ets grow initially very fast before settling into an os- 
cillatory quasi-static state which changes slowly on a 
secular time scale. This slow increase of the eccentric- 
ities on the longer time scale is caused by the growing 
gravitational forces between the planets, due to the 
decreasing radial distance of the two planets on their 
inward migration process. 

c) The ellipses/periastrae of the planets rotate at a con- 
stant, retrograde angular speed w. Coupled together 
by the resonance, the apsidal precession rate w for 
both planets is identical, which can be inferred from 
the parallel lines in the bottom panel of Fig. 2. The 



orientation of the orbits is phase-locked with a con- 
stant separation Aw = tu 2 — vj\ . The rotation period 
of the ellipses (apsidal lines) is slightly longer than the 
oscillation period of the eccentricities. 
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Fig. 5. Eccentricity ratio e2/ei of the outer and inner planet, 
versus periastron difference Aw for model C. The data points 
are spaced equally in time with a distance of approximately 
St = 23 yrs. Shown is a section of the evolution of model C, 
from 26,000 to 55,000 yrs, which covers about 7 libration pe- 
riods. 
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Time (Years) 

Fig. 6. The semi-major axis, eccentricity and position angle 
of the orbital periastron versus time for Model D. The only 
parameter different from the first model C is the lower viscos- 
ity. Here, the outer planet passes through the 3:1 resonance 
and captures the inner planet finally into a 2:1 configuration. 
The dotted reference line (labeled 2:1) marks the location of 
the 2:1 resonance with respect to the inner planet. Upon reso- 
nant capture, the eccentricities grow and the two orbits librate 
retrograde with a fixed relative orientation of Azu — 0°. 



The capture into resonance and the subsequent libra- 
tion of the orbits in model C is illustrated further in Fig. 4. 
As suggested in Fig. 2 (bottom panel) the periastrae be- 
gin to align upon capture in the 3:1 resonance. Initially, 
during the phase when the eccentricities are still rising 
(between 5 and 20 thousand yrs), the difference of the pe- 
riastrae settles intermediately to Aw w 180° (see Fig. 4, 
left panel). Then, upon saturation after about 20,000 yrs, 
the system re-adjusts and eventually establishes itself at 
Aw « 107°, with a libration amplitude of about 7°. The 
right panel of Fig. 4 shows the time evolution of the reso- 
nant angle 0i = 3A2 — Ai — 2w\. Here Ai and A2 denote 
the mean longitudes of the inner and outer planet, respec- 
tively. Initially, 0i settles to 180° as well, and re-adjusts 
after 20,000 yrs to —145° with a libration amplitude of 
±15°. This behavior of an initially symmetric alignment 
of Aw and Qi at about 180° followed by a later change 
to Aw w ±110° and 81 w ±145° is characteristic for all 



our models which show captures into 3:1 resonance, inde- 
pendent of the physical parameters (see Tab. 4). 

This behavior can be understood by an analysis of 
the interaction Hamiltonian for resonant systems (Beauge 
et al. 2002). By minimizing the interaction energy, the 
equilibrium values for Aw and 61 can be obtained as a 
function of the mass and eccentricity of the two planets. 
As shown by Beauge et al. (2002), when the eccentricity of 
the inner planet is small (ei <, .12) the equilibrium values 
of both resonant angles are exactly 180°. For higher eccen- 
tricity though, the equilibrium values of Aw and 0i shift 
to 115° and 210°. In our numerical simulations we find ex- 
actly this behavior. Initially, upon entering the resonant 
configuration the eccentricities are small and the two an- 
gles both adjust to 180°. Later they readjust to new values 
as the eccentricities rise above the critical value. 

In the subsequent longterm evolution after 20,000 yrs, 
the system settles into a quasi-equilibrium situation where 
the eccentricities oscillate with a period of about 3750 yrs. 
In the 62/ ei versus Aw diagram (Fig. 5), this phenomenon 
is demonstrated by the circular distribution of data points 
around the equilibrium values. 

4.1.3. 2:1 Resonance: Model D 

In comparison to the previous model C the only difference 
in this model D is the value of the viscosity coefficient. 
Three times less viscosity (a = 3.3- 10~ 3 ) results in a little 
bit slower migration speed. For this model the evolution 
did not end up in the 3:1, but rather the 2:1, resonance. In 
Figure 6, the evolution of a, e and w is displayed. The ec- 
centricities show a small 'kink' as the outer planet reaches 
the 3:1 resonance, but migration continues past the reso- 
nant location. Later, at about 26,000 yrs, capture occurs 
in the 2:1 resonance, leading to perfectly aligned orbits. 
Both resonant angles, Aw and Oi, are nearly zero with a 
very small libration amplitude. The eife\ versus Aw plot 
(Fig. 8) shows the small variations in the eccentricities and 
the small libration amplitude (cf. Fig. 5). Using the same 
model parameters, we ran a nearly identical simulation in 
which the initial density was disturbed randomly by 1%. 
With just this small change in the initial conditions, cap- 
ture into 3:1 resonance was successful (see also Tab. 4 and 
Sect. 4.3.3). 

4.2. Damped N-Body Calculations 

As a test of the damped N-body model described above 
(Sect. 3.2), we calculate the evolution of all models C- 
G using the 3-body method and compare the results to 
the full hydrodynamic evolutions. As outlined above, the 
prescribed damping formula (Eq. 6) is used to directly 
alter the semi-major axis a and eccentricity e of the outer 
planet only. 

All the results of the damped N-body models are sum- 
marized in Tab. 3. The values of the damping constant 
K cannot be determined too precisely, as there is some 
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Fig. 7. Results for model D. Left: The difference in the direction of the periastrae, Acc = VJ2 — wi (in degrees) for the two 
planets versus time (in thousand years). Right: The resonant angle for the 2:1 resonance Oi = 2A 2 — Ai — zui (in degrees), as 
a function of time. 



Table 3. Fit parameters of the N-body computations, obtained 
through comparison with the full hydro simulations for mod- 
els C-G. Listed are the model name, the initial damping time 
scale ro, the slow down of the damping (3 (see Eq. 6), and the 
eccentricity damping K (see Eq. 8). 



Model 


to [yrs] 


P 


K 


C 


17,500 


3.2 


1.5 


D 


27,000 


2.5 


1.5 


E 


38,000 


2.5 


1.5 


F 


19,500 


1.5 


1.5 


G 


33,500 


1.0 


2.5 



dependence of the magnitude on the initial eccentricities. 
Those were chosen to be about 0.01 to 0.02 for all models. 

As an example, we display the evolution of model G, 
in which H/r = 0.05. To test the damping, we first use 
Eq. (6) with a constant damping r = to, i.e. (3 = 0. This 
refers to the lower (dark gray) curve in Fig. (9), where 
we chose t = 33, 500 yrs. As can be seen in the figure, 
a constant damping, even when it has the correct initial 
slope, in the longterm yields too fast a migration of the 
outer planet. A time-dependent damping with (3 = 1.0 
(light gray curve), leads to a much better fit with the 
hydrodynamic results. Additionally, we tested our method 
on published results for migrating single planets (Nelson 
et al. 2000), and find good agreement for suitable To and 
0. Pure exponential fits for 02 (t) with (3 = generally do 
not give satisfactory results. 

This decrease in the damping rate as a function of time 
is a result of the reduction of mass in the disk, mainly due 
to the accretion of matter onto the outer planet. The mass 
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Fig. 8. Eccentricity ratio ez/ei of the outer and inner planet, 
versus periastron difference Am for model D. The data points 
are spaced equally in time with a distance of approximately 
St — 22 yrs. Shown is a section of the evolution of model D, 
from 43,000 to about 74,000 yrs. 

flow across the gap is small, and the accretion rate onto 
the inner planet is substantially lower (see Kley 2000). 
However, if no gas is allowed to accrete onto the outer 
planet, a more substantial amount of gas may flow oc- 
cur across the gap, leading to a smaller migration rate, 
or even outward migration (Masset & Snellgrove 2001; 
Masset 2002; Masset & Papaloizou 2003), because the an- 
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Fig. 9. The evolution of the semi-major axis of the outer planet 
a.2, comparing a hydrodynamic model with two N-body calcu- 
lations. The black curve (labeled Hydro) is the result of the 
full hydro model G. The lower dark gray curve is the damped 
N-body model using a constant damping r = to in Eq. 6. The 
upper lighter gray curve (on top of the black one) is for time- 
dependent damping with a non-zero value of (3 = 1.0, resulting 
in a much better match with the hydro model. 



gular momentum lost by the gap crossing material will be 
gained by the planet. 

Besides a, we also compared the evolution of the ec- 
centricity e between the two approaches. In Fig. (10) we 
display the eccentricity evolution of the full hydro and 
the damped 3-body case for model G. Despite some dif- 
ferences which we attribute to the unknown eccentricity 
damping mechanism, the overall agreement is reasonable. 
For a given semi-major axis damping rate, the final values 
obtained for e\ and e-i at larger times depend on the initial 
values for the eccentricities and the amount of eccentricity 
damping. In this case we used an initial e(t Q ) = 0.01 for 
both planets and an eccentricity damping factor K — 2.5, 
i.e. a slightly shorter damping time scale for eccentricity as 
for semi-major axis. For all models we find that the eccen- 
tricity damping rate is of the same order as the semi-major 
axis damping, i.e. K = 0(1). This finding is in contrast 
to Lee & Peale (2002a) who determined a much shorter 
eccentricity damping time, based on models for GJ 876. 

There are several possible reasons for this difference. 
Since the eccentricity damping of a planet is caused by 
material in the co-orbital region close to the planet, the 
treatment of this region in the models may have some 
effect on the results. In particular, the mass accretion onto 
the planet, the smoothing of the gravitational potential, 
and the numerical resolution may each play a significant 
role here. However, the simulations by Lee & Peale (2002a) 
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Fig. 10. The evolution of the eccentricities of the inner planet 
(upper curves) and the outer planet (lower curves) for the full 
hydro model G (black curves covering the whole time range) 
and the damped N-body model (light gray), using a damping 
constant of K = 2.5. 
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Fig. 11. The time evolution of the outer planet's semi-major 
axis for three models (C-E) with different viscosities. The cor- 
responding a values are given in the legend 



are not based on clear physical model of the damping, but 
rather use an ad hoc prescription. Fitting to the observed 
case of GJ 876 yields a high value of K. On the other 
hand, it will be difficult to model the system HD 82943 
with K — 100 because of the high observed eccentricities. 
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4.3. Dependence on the Physical Parameters 

After describing the major effects of resonant capture and 
evolution we focus now on the dependence on the physical 
parameters. An overview of the results for all models is 
given in Tab. 4. 

Table 4. Results of the full hydrodynamical computations for 
the longterm evolution. Listed are the model name, the reso- 
nance in which the system is captured, the eccentricity ratio 
e2/ei, the separation of the periastrae Azu, the resonant angle 
©i, and the speed of the apsidal precession zu. The models la- 
beled with (V) are still evolving with time. For model D two 



cases (cau 


ght in 


different 


resonances) are 


presented. 


Model 


Res. 


e2/ei 


Azu 


©i 


zb 








[deg] 


[deg] 


[rad/yr] 


A 


3:1 


1.0 


-110 


-140 


- 0.0015 


B 


3:1 


0.7 


-120 


-150 


- 0.0015 


C 


3:1 


0.73 


-107 


-145 


- 0.0014 


Dl 


2:1 


0.35 








- 0.0033 


D2 


3:1 


0.70 


+110 


147 


- 0.0015 (V) 


E 


3:1 


0.73 


-107 


-140 


- 0.0010 


F 


5:2 


1.1 


+180 





- 0.0033 


G 


3:1 


0.65 


+110 


+150 


- 0.0022 (V) 


X 


2:1 


0.24 








- 0.0021 
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Fig. 12. The evolution of the semi-major axis of the outer 
planet ai for two models (C, G) with different vertical disk 
scale heights, as indicated in the legend. 



4.3.1. Viscosity 

As the disk viscosity, parameterized here through the stan- 
dard a-parameter, determines the overall evolution of the 
disk, it is to be expected that its magnitude influences 
the longterm evolution of the planetary system as well. 
Starting from the standard model C (a = 10~ 2 ) we per- 
formed additional runs using different values for the vis- 
cosity parameter a; model D has a — 3.3- 10~ 3 and model 
E has a = 10~ 3 . The semi-major axis evolution for these 
models is shown in Fig. 11. For a single planet, the migra- 
tion rate depends on the value of the viscosity (Nelson 
et al. 2000). Hence, the initial migration rates (during 
the first 10,000 yrs) are reduced for smaller v. The ini- 
tial migration speeds for each model, as seen in Fig. 11, 
are quantified as tq in Tab. 3. 

Upon capture of the inner planet, the migration rate 
of the outer slows down. Contrary to the expectation that 
smaller viscosities yield slower migration, we find that for 
these resonantly driven double-planet systems, the speed 
of migration is eventually faster for smaller viscosity coef- 
ficients. When the viscosity is lower, the migration slows 
down less rapidly (this would correspond to a lower (3- 
valuc in the N-body method). This time-dependence of 
the damping is linked to the assumed mass accretion onto 
the planet. The accretion, which is larger for higher vis- 
cosity, lowers the mass of the disk with time and reduces 
the effective torques, such that migration drops off more 
rapidly for higher viscosity. The ratio of semi-major axis 
damping to eccentricity damping (K), however, does not 
appear to depend on the value of the viscosity. In all three 



cases we find capture into the 3:1 resonance (see Tab. 4), 
although this result is somewhat indeterminate. For the 
value of a = 3.3 • 10~ 3 , we also find a case in 2:1. 

4.3.2. Temperature 

Similarly to the viscosity, it may be expected that the 
scale height (temperature) of the disk will also influence 
the migration process of planets. Results for single plan- 
ets (Kley 1999) have shown that for larger H/r the gap is 
less cleared because the larger pressure gradient 'pushes' 
material into the gap. The higher density near the gap 
edges leads to a faster migration speed. This effect is in- 
deed clearly seen in the early evolution of the planetary 
systems C & G, where the decline in the semi-major axis 
of the outer planet a 2 is faster for higher temperatures 
(see Fig. 12). 

Upon capture however, the evolution behaves similarly 
to the case of varying viscosity. Those systems with higher 
H/r now have the lowest migration speed. Similarly to the 
models with lower viscosity, smaller H/r leads to lower ac- 
cretion, higher disk mass, and hence more material push- 
ing on the planets. As the migration slows down less for 
lower H/r, the overall evolution time scale becomes so 
long that it takes much more than 10 5 yrs for the sys- 
tem to reach an equilibrium state. The intermediate model 
with H/R = 0.075 is not displayed as is went into a 5:2 
resonance. 



12 



Kley, Peitz, and Bryden: Planetary Systems in Resonance 



200 
150 
100 



„ 50 
s 



-50 
-100 
-150 
-200 



_-0--0^-O ^ ^ 




0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 
m p [ M Jup] 



Fig. 13. Results of a sequence of damped N-body simulations. 
Plotted is the difference of the periastrae Azu — zu2 — rui of the 
two planets after capture into resonance versus planet mass, 
where m p = m2 = mi. The diamonds indicate capture in 2:1 
resonance, while the plus signs are for 3:1 resonance. The other 
parameters are fixed, as described in the tex. The horizontal 
lines indicate values of 0, ±100° and ±180°. 

4.3.3. Capture Probabilities 

The main results of the full hydrodynamic calculations, 
as summarized in Tab. 4, show that for the same physi- 
cal setup, capture in different resonances may occur. As a 
three body system is intrinsically chaotic, this indetermi- 
nate behavior may be expected. 

Nevertheless, by running a whole sequence of fast 
damped N-body models, we can investigate what condi- 
tions determine the principle final outcome. The standard 
setup consists of two planets of lMj up each, placed ini- 
tially at 4 and 12 AU from the central 1M© star. The ini- 
tial eccentricities are 0.02. As in all previous models, only 
the orbit of the outer planet is damped, using in these 
cases a damping time scale To = 20, 000 yrs. We fix the 
damping constants to (3 = 1.0 and K = 1.0. Starting from 
this standard case, we vary the damping time scale To, the 
initial eccentricity e 2 , and mass of the outer planet m 2 , 
while keeping always identical planet masses, mi = m 2 . 

In the standard case, the planets are caught in a 2:1 
resonance. Varying the time scale To from 10,000 to 50,000 
yrs, there is no capture into higher resonances. Upon vari- 
ation of eccentricity e 2 from 0.0 to 0.5, we find that for 
e 2 < 0.25 capture occurs always into 2:1, while for larger 
e 2 higher order resonances (primarily 3:1) are possible, 
however with no definite outcome. The influence of the 
planet mass is illustrated in Fig. 13, where the resonance 
type (2:1 diamonds, 3:1 plus signs) is shown as a function 
of planet mass (with mi = m 2 ). For small planets with 
m p < 1.7Mj up , capture occurs robustly into the 2:1 reso- 



nance. Higher resonances are possible only for m p larger 
than this value. However, due to the chaotic nature of 
the problem the exact outcome for a particular m p is not 
predictable. This is in agreement with the hydrodynamic 
simulations where we also find capture in different reso- 
nances just by perturbing the initial density slightly (com- 
pare models Dl, D2). 



4.3.4. The Resonant Angles 

In the majority of hydrodynamic models, the planets catch 
each other in a 3:1 resonance. As demonstrated above, 
this is mainly a consequence of the large masses chosen 
for the planets. For all models the resonant angles settle 
to \Azu\ « 110° and |©i| w 145°, with some scatter of 
about ±10°. For the damped N-body models we checked 
the values of the resonant angles as well, finding that all 
2:1 resonances settle into the complete symmetric config- 
uration Aw = 61 = (Fig. 13). For the 3:1 resonances 
we see anti-symmetric configurations with anti-aligned pe- 
riastrae, I Aw I = 180° and 81 = for masses around 
m p — 2Mj up , (when the first 3:1 cases begin to occur), 
while for all larger masses we find preferentially the previ- 
ous non-symmetric configurations. This behavior is again 
an indication of a bifurcation in the stability properties 
of resonant systems, as claimed by Beauge et al. (2002). 
With additional damped N-body calculations we also find 
that systems entering a 5:2 resonance (not shown) ex- 
hibit a typical anti-symmetric behavior \Aw\ = 180° and 
6i = 0°. 

5. Summary and Conclusion 

We have performed full hydrodynamical calculations sim- 
ulating the joint evolution of a pair of protoplanets to- 
gether with the surrounding protoplanctary disk from 
which they originally formed. The focus lies on massive 
planets in the range of a few Jupiter masses. For the disk 
evolution we solve the Navier-Stokes equations, and the 
motion of the planets is followed using a 4th order Runge- 
Kutta scheme, which includes their mutual interactions as 
well as the star and disk's gravitational fields. These re- 
sults were compared to simplified (damped) N-body com- 
putations, where the gravitational influence of the disk is 
modeled through analytic damping terms applied to the 
semi-major axis and eccentricity. 

We find that both methods yield comparable results, 
if the damping constants in the simplified models are ad- 
justed properly. The mass reduction of the disk with time, 
due for example to mass accretion onto the planet, or pos- 
sible mass flow across the outer planet's gap can be mod- 
eled satisfactorily through a damping time scale, which 
depends linearly on time. The eccentricity damping was 
always chosen to be a constant multiple K of the semi- 
major axis damping. In this case we find that K must be 
of order unity to match the hydrodynamic models. 

However, fitting N-body models to the observed pa- 
rameters of GJ 876 requires a high e-damping with typ- 
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ically K = 100 (see also Lee & Peale 2002a), relatively 
independent of the functional behavior of a<i (t). Reasons 
for this discrepancy may lie in the simplified hydrodynam- 
ical model, which uses a fixed equation of state, a simple 
treatment of the planetary structure, and only an approx- 
imate model of the torques acting on the planet. Also, 
eccentricity damping is dominated by material close to 
the planet; the insufficient numerical grid resolution near 
the planet may smear out the damping forces. In addi- 
tion, the accretion process of matter onto the planet is 
reducing the mass in the co-orbital region which lowers 
the eccentricity damping. The simplified assumption of a 
constant value of K needs to be checked. More detailed 
hydrodynamical models may help to resolve this discrep- 
ancy in the future. An alternative explanation for the low 
eccentricities in GJ 876 compared to our hydrodynamic 
simulations is that further evolution of the eccentricities 
occurs in the system after planet and gap formation. The 
planet eccentricities may be further modified as the disk 
dissipates and its resulting eccentricity forcing gradually 
declines. On the other hand, the assumption of a constant 
value of K in the N-body models in the computations by 
Lee & Peale (2002a) is also not based on any detailed hy- 
drodynamic model but rather assumed ab initio. In more 
general models, this will have to be relaxed. 

The case HD 82943 is also not easy to model as the ec- 
centricities for both planets are very large, which turned 
out to be very difficult to capture with N-body models, 
even with very low damping. The problem here lies in the 
stability of the resulting system. All test computations 
with constant values of K eventually led to unstable sys- 
tems. Compared to GJ 876, the eccentricity damping for 
HD 82943 must be orders of magnitude less, if otherwise 
similar physical parameters are used. In order to explain 
the high eccentricities, the inclusion of an additional com- 
panion may be necessary. 

Despite of the difficulty of the models to obtain the 
observed eccentricities, there are nevertheless several fea- 
tures of the observed 2:1 planets which are captured cor- 
rectly by our simulations: i) The larger mass of the outer 
planet, it) the higher eccentricity of the inner planet, and 
Hi) the periastrae separation of Am = 0°. These are ro- 
bust predictions of the hydrodynamic models. 

For 3:1 resonances, anti-symmetric (Aw = 180°) and 
non-symmetric final configurations are obtained. In the 
non-symmetric case we found over a range of models a 
value of | Aw ~ 110°|, which is supported by stability anal- 
ysis (Beauge et al. 2002). In 55 Cnc, the only observed 3:1 
case, there are other planets present in the system, which 
makes an interpretation using just this simple treatment 
questionable. 
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